This question involves the use of simple linear regression on the Auto data set.
library(data.table)
library(tidyverse)
auto <- fread("../Assignment 1/Auto.csv")
auto$horsepower <- as.numeric(auto$horsepower)
Warning: NAs introduced by coercion
Use the
lm()function to perform a simple linear regression withmpgas the response andhorsepoweras the predictor. Use thesummary()function to print the results. Comment on the output. For example:
lm <- lm(mpg ~ horsepower, data = auto)
lm
Call:
lm(formula = mpg ~ horsepower, data = auto)
Coefficients:
(Intercept) horsepower
39.9359 -0.1578
summary(lm)
Call:
lm(formula = mpg ~ horsepower, data = auto)
Residuals:
Min 1Q Median 3Q Max
-13.5710 -3.2592 -0.3435 2.7630 16.9240
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.935861 0.717499 55.66 <2e-16 ***
horsepower -0.157845 0.006446 -24.49 <2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.906 on 390 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.6059, Adjusted R-squared: 0.6049
F-statistic: 599.7 on 1 and 390 DF, p-value: < 2.2e-16
Is there a relationship between the predictor and the re- sponse?
Yes. After horsepower is converted to numeric, there is a relationship between horsepower and mpg.
How strong is the relationship between the predictor and the response?
Very strong. The p-value is tiny.
Is the relationship between the predictor and the response positive or negative?
Negative. As horsepower increases, mpg decreases.
What is the predicted
mpgassociated with ahorsepowerof 98? What are the associated 95% confidence and prediction intervals?
predict(lm, auto)[98]
98
23.36216
# confidence interval
predict(lm, auto, interval = "confidence")[98,]
fit lwr upr
23.36216 22.87497 23.84936
# prediction interval
predict(lm, auto, interval = "prediction")[98,]
fit lwr upr
23.36216 13.70483 33.01950
Plot the response and the predictor. Use the
abline()function to display the least squares regression line.
plot(auto$horsepower, auto$mpg)
abline(lm)
Use the
plot()function to produce diagnostic plots of the least squares regression fit. Comment on any problems you see with the fit.
plot(lm)
From the textbook, Figure 3.9: “Plots of residuals versus predicted (or fitted) values for the Auto data set. In each plot, the red line is a smooth fit to the residuals, intended to make it easier to identify a trend. Left: A linear regression of mpg on horsepower. A strong pattern in the residuals indicates non-linearity in the data. Right: A linear regression of mpg on horsepower and horsepower. There is little pattern in the residuals.”
This question involves the use of multiple linear regression on the
Autodata set.
Produce a scatterplot matrix which includes all of the variables in the data set.
auto$name <- as.factor(auto$name)
pairs(auto)
Compute the matrix of correlations between the variables using the function
cor(). You will need to exclude thenamevariable, which is qualitative.
cor(auto[, -9])
mpg cylinders displacement horsepower
mpg 1.0000000 -0.7762599 -0.8044430 NA
cylinders -0.7762599 1.0000000 0.9509199 NA
displacement -0.8044430 0.9509199 1.0000000 NA
horsepower NA NA NA 1
weight -0.8317389 0.8970169 0.9331044 NA
acceleration 0.4222974 -0.5040606 -0.5441618 NA
year 0.5814695 -0.3467172 -0.3698041 NA
origin 0.5636979 -0.5649716 -0.6106643 NA
weight acceleration year origin
mpg -0.8317389 0.4222974 0.5814695 0.5636979
cylinders 0.8970169 -0.5040606 -0.3467172 -0.5649716
displacement 0.9331044 -0.5441618 -0.3698041 -0.6106643
horsepower NA NA NA NA
weight 1.0000000 -0.4195023 -0.3079004 -0.5812652
acceleration -0.4195023 1.0000000 0.2829009 0.2100836
year -0.3079004 0.2829009 1.0000000 0.1843141
origin -0.5812652 0.2100836 0.1843141 1.0000000
Use the
lm()function to perform a multiple linear regression withmpgas the response and all other variables exceptnameas the predictors. Use thesummary()function to print the results. Comment on the output. For instance:
lm <- lm(mpg ~ cylinders + displacement + horsepower + weight + acceleration + year + origin, auto)
lm
Call:
lm(formula = mpg ~ cylinders + displacement + horsepower + weight +
acceleration + year + origin, data = auto)
Coefficients:
(Intercept) cylinders displacement horsepower
-17.218435 -0.493376 0.019896 -0.016951
weight acceleration year origin
-0.006474 0.080576 0.750773 1.426140
summary(lm)
Call:
lm(formula = mpg ~ cylinders + displacement + horsepower + weight +
acceleration + year + origin, data = auto)
Residuals:
Min 1Q Median 3Q Max
-9.5903 -2.1565 -0.1169 1.8690 13.0604
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.218435 4.644294 -3.707 0.00024 ***
cylinders -0.493376 0.323282 -1.526 0.12780
displacement 0.019896 0.007515 2.647 0.00844 **
horsepower -0.016951 0.013787 -1.230 0.21963
weight -0.006474 0.000652 -9.929 < 2e-16 ***
acceleration 0.080576 0.098845 0.815 0.41548
year 0.750773 0.050973 14.729 < 2e-16 ***
origin 1.426141 0.278136 5.127 4.67e-07 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.328 on 384 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.8215, Adjusted R-squared: 0.8182
F-statistic: 252.4 on 7 and 384 DF, p-value: < 2.2e-16
Is there a relationship between the predictors and the re- sponse?
There is a relationship between mpg and displacement, weight, year, and origin.
Which predictors appear to have a statistically significant relationship to the response?
There is a relationship between mpg and displacement, weight, year, and origin.
What does the coefficient for the
yearvariable suggest?
This suggests that the mpg increases alongside the year. (Cars become more efficient over time).
Use the
plot()function to produce diagnostic plots of the linear regression fit. Comment on any problems you see with the fit. Do the residual plots suggest any unusually large outliers? Does the leverage plot identify any observations with unusually high leverage?
plot(lm)
There are a decent amount of outliers among the dataset. There are a few observations with strangely high leverage.
Use the
*and:symbols to fit linear regression models with interaction effects. Do any interactions appear to be statistically significant?
summary(lm(mpg ~ cylinders * displacement + horsepower:weight, auto))
Call:
lm(formula = mpg ~ cylinders * displacement + horsepower:weight,
data = auto)
Residuals:
Min 1Q Median 3Q Max
-15.8087 -2.2841 -0.5424 2.1259 17.5877
Coefficients:
Estimate Std. Error t value
(Intercept) 5.524e+01 2.384e+00 23.167
cylinders -3.831e+00 5.339e-01 -7.176
displacement -1.388e-01 1.511e-02 -9.184
cylinders:displacement 1.923e-02 2.171e-03 8.857
horsepower:weight -2.231e-05 2.961e-06 -7.533
Pr(>|t|)
(Intercept) < 2e-16 ***
cylinders 3.68e-12 ***
displacement < 2e-16 ***
cylinders:displacement < 2e-16 ***
horsepower:weight 3.54e-13 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.164 on 387 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.7182, Adjusted R-squared: 0.7153
F-statistic: 246.6 on 4 and 387 DF, p-value: < 2.2e-16
I guessed and it turns out that my suspicions were correct. cylinders:diplacement and horsepower:weight were significant.
Try a few different transformations of the variables, such as \(log(X)\), \(\sqrt{X}\), \(X^2\). Comment on your findings.
summary(lm(log(mpg) ~ cylinders^2, auto))
Call:
lm(formula = log(mpg) ~ cylinders^2, data = auto)
Residuals:
Min 1Q Median 3Q Max
-0.61694 -0.12149 -0.00995 0.12358 0.62573
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.002754 0.032379 123.62 <2e-16 ***
cylinders -0.165149 0.005664 -29.16 <2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1918 on 395 degrees of freedom
Multiple R-squared: 0.6828, Adjusted R-squared: 0.682
F-statistic: 850.3 on 1 and 395 DF, p-value: < 2.2e-16
summary(lm(log(mpg) ~ cylinders^2 + log(horsepower)^2:sqrt(weight), auto))
Call:
lm(formula = log(mpg) ~ cylinders^2 + log(horsepower)^2:sqrt(weight),
data = auto)
Residuals:
Min 1Q Median 3Q Max
-0.47818 -0.10004 -0.00615 0.10030 0.51869
Coefficients:
Estimate Std. Error t value
(Intercept) 4.4881118 0.0421172 106.563
cylinders -0.0211598 0.0107477 -1.969
log(horsepower):sqrt(weight) -0.0050924 0.0003451 -14.756
Pr(>|t|)
(Intercept) <2e-16 ***
cylinders 0.0497 *
log(horsepower):sqrt(weight) <2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1535 on 389 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.7972, Adjusted R-squared: 0.7961
F-statistic: 764.4 on 2 and 389 DF, p-value: < 2.2e-16
I’m a little surprised that the regressions didn’t break as easily as I thought they would with some basic transformations.
This question should be answered using the
Carseatsdata set.
carseats <- ISLR2::Carseats
Fit a multiple regression model to predict
SalesusingPrice,Urban, andUS.
lm <- lm(Sales ~ Price + Urban + US, carseats)
lm
Call:
lm(formula = Sales ~ Price + Urban + US, data = carseats)
Coefficients:
(Intercept) Price UrbanYes USYes
13.04347 -0.05446 -0.02192 1.20057
summary(lm)
Call:
lm(formula = Sales ~ Price + Urban + US, data = carseats)
Residuals:
Min 1Q Median 3Q Max
-6.9206 -1.6220 -0.0564 1.5786 7.0581
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.043469 0.651012 20.036 < 2e-16 ***
Price -0.054459 0.005242 -10.389 < 2e-16 ***
UrbanYes -0.021916 0.271650 -0.081 0.936
USYes 1.200573 0.259042 4.635 4.86e-06 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.472 on 396 degrees of freedom
Multiple R-squared: 0.2393, Adjusted R-squared: 0.2335
F-statistic: 41.52 on 3 and 396 DF, p-value: < 2.2e-16
Provide an interpretation of each coefficient in the model. Be careful—some of the variables in the model are qualitative!
As Price increases by 1, Sales decrease by \(.054459\) When Urban is Yes, Sales decrease by \(.021916\) (but this isn’t significant) When US is Yes, Sales increases by \(1.021916\)
Write out the model in equation form, being careful to handle the qualitative variables properly.
\(\textrm{Sales}=13.04-.05\textrm{Price}-(.02\textrm{Urban}_{\textrm{Yes}})+(1.20\textrm{US}_{\textrm{Yes}})\)
For which of the predictors can you reject the null hypothesis \(H0:β_j=0\) ?
Price and USYes are significant enough with a p-value less than 0.05 so we are 95% confident that they are significant.
On the basis of your response to the previous question, fit a smaller model that only uses the predictors for which there is evidence of association with the outcome.
lm <- lm(Sales ~ Price + US, carseats)
lm
Call:
lm(formula = Sales ~ Price + US, data = carseats)
Coefficients:
(Intercept) Price USYes
13.03079 -0.05448 1.19964
summary(lm)
Call:
lm(formula = Sales ~ Price + US, data = carseats)
Residuals:
Min 1Q Median 3Q Max
-6.9269 -1.6286 -0.0574 1.5766 7.0515
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.03079 0.63098 20.652 < 2e-16 ***
Price -0.05448 0.00523 -10.416 < 2e-16 ***
USYes 1.19964 0.25846 4.641 4.71e-06 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.469 on 397 degrees of freedom
Multiple R-squared: 0.2393, Adjusted R-squared: 0.2354
F-statistic: 62.43 on 2 and 397 DF, p-value: < 2.2e-16
How well do the models in (a) and (e) fit the data?
The model in (e) fits much better than the one in (a) because it only has significant variables.
Using the model from (e), obtain \(95\%\) confidence intervals for the coefficient(s).
confint(lm)
2.5 % 97.5 %
(Intercept) 11.79032020 14.27126531
Price -0.06475984 -0.04419543
USYes 0.69151957 1.70776632
Is there evidence of outliers or high leverage observations in the model from (e)?
par(mfrow = c(2,2))
plot(lm)
There are a few points that look to be outside the normal range. Points 51, 69, and 377 seem to be outliers and points 26, 50, and 368 seem to have high leverage. I’m not worried about them because there are only three in a good-sized dataset.
From the data available at http://www.stat.ufl.edu/~winner/datasets.html, obtain the data on “Fibre Diameters and Breaking Strenghs for Nextel 610 Fibres.” (please note that there is a typo on the website. It should be Strength and not Strengh) According to the description available there, the expectation is that the log of breaking strength of the fibre should be negatively and linearly related to diameter. (Note log here means natural log if not specified.)
dt <- fread("fibre.csv")
colnames(dt) <- c("stren", "diam")
head(dt)
Produce a scatter plot of breaking strength against diameter.
plot(dt$diam, dt$stren)
Produce a scatter plot of the log of breaking strength against diameter.
dt$logStren <- log(dt$stren)
plot(dt$logStren, dt$diam)
Produce a scatter plot of the log of breaking strength against the log of diameter.
dt$logDiam <- log(dt$diam)
plot(dt$logStren, dt$logDiam)
Regress breaking strength on diameter.
lm <- lm(stren ~ diam, data = dt)
lm
Call:
lm(formula = stren ~ diam, data = dt)
Coefficients:
(Intercept) diam
5312.8 -197.7
summary(lm)
Call:
lm(formula = stren ~ diam, data = dt)
Residuals:
Min 1Q Median 3Q Max
-603.46 -246.10 -6.76 223.64 1047.22
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5312.78 1040.41 5.106 5.61e-06 ***
diam -197.73 91.91 -2.151 0.0365 *
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 335.7 on 48 degrees of freedom
Multiple R-squared: 0.08794, Adjusted R-squared: 0.06894
F-statistic: 4.628 on 1 and 48 DF, p-value: 0.03651
Regress the log of breaking strength on diameter.
lm <- lm(logStren ~ diam, data = dt)
lm
Call:
lm(formula = logStren ~ diam, data = dt)
Coefficients:
(Intercept) diam
8.76091 -0.06504
summary(lm)
Call:
lm(formula = logStren ~ diam, data = dt)
Residuals:
Min 1Q Median 3Q Max
-0.211992 -0.075119 0.003561 0.076367 0.293763
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.76091 0.33575 26.094 <2e-16 ***
diam -0.06504 0.02966 -2.193 0.0332 *
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1083 on 48 degrees of freedom
Multiple R-squared: 0.09105, Adjusted R-squared: 0.07211
F-statistic: 4.808 on 1 and 48 DF, p-value: 0.0332
Regress the log of breaking strength on the log of diameter.
lm <- lm(logStren ~ logDiam, data = dt)
lm
Call:
lm(formula = logStren ~ logDiam, data = dt)
Coefficients:
(Intercept) logDiam
9.8328 -0.7454
summary(lm)
Call:
lm(formula = logStren ~ logDiam, data = dt)
Residuals:
Min 1Q Median 3Q Max
-0.21123 -0.07610 0.00396 0.07673 0.29398
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.8328 0.8209 11.978 4.99e-16 ***
logDiam -0.7454 0.3385 -2.202 0.0325 *
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1083 on 48 degrees of freedom
Multiple R-squared: 0.09174, Adjusted R-squared: 0.07282
F-statistic: 4.848 on 1 and 48 DF, p-value: 0.03251
From the data available at http://www.stat.ufl.edu/~winner/datasets.html, obtain the data on “Variables associated with Permeability and Porosity of Rocks”
dt <- fread("rocks.csv")
colnames(dt) <- c("type", "density", "porosity", "logPermeability", "residue", "carbonate", "aGrain", "aSD", "bGrain", "bSD", "calcite", "dolomite")
head(dt)
Fit a multiple regression model to predict porosity. Please provide a clean model with only significant variables. Interpret your model results and diagnostics.
lm <- lm(porosity ~ ., data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
# This prunes the regression until only significant variables remain
while(length(coefs) != sum(coefs)) {
form <- paste("porosity ~", paste(names(coefs[coefs]), "+", collapse = " "), collapse = " ")
form <- substr(form, 1, nchar(form) - 2)
form
lm <- lm(formula = form, data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
}
lm
Call:
lm(formula = form, data = dt)
Coefficients:
(Intercept) logPermeability
6.240 2.666
lmSumm
Call:
lm(formula = form, data = dt)
Residuals:
Min 1Q Median 3Q Max
-2.1522 -1.0096 -0.2271 0.3610 5.8290
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.240 2.046 3.050 0.00485 **
logPermeability 2.666 1.222 2.181 0.03740 *
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.689 on 29 degrees of freedom
Multiple R-squared: 0.141, Adjusted R-squared: 0.1113
F-statistic: 4.759 on 1 and 29 DF, p-value: 0.0374
par(mfrow = c(2,2))
plot(lm)
Fit a multiple regression model to predict log(permeability). Please provide a clean model with only significant variables. Interpret your model results and diagnostics.
lm <- lm(logPermeability ~ ., data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
# This prunes the regression until only significant variables remain
while(length(coefs) != sum(coefs)) {
form <- paste("logPermeability ~", paste(names(coefs[coefs]), "+", collapse = " "), collapse = " ")
form <- substr(form, 1, nchar(form) - 2)
form
lm <- lm(formula = form, data = dt)
lmSumm <- summary(lm)
coefs <- coef(lmSumm)[-1,4] <= .05
}
lm
Call:
lm(formula = form, data = dt)
Coefficients:
(Intercept) porosity aGrain bGrain
-0.448085 0.063658 0.897464 -0.994242
calcite
0.003091
lmSumm
Call:
lm(formula = form, data = dt)
Residuals:
Min 1Q Median 3Q Max
-0.58844 -0.06918 0.02635 0.09975 0.40127
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.448085 0.499277 -0.897 0.37771
porosity 0.063658 0.023650 2.692 0.01227 *
aGrain 0.897464 0.430584 2.084 0.04710 *
bGrain -0.994242 0.431684 -2.303 0.02952 *
calcite 0.003091 0.001066 2.899 0.00751 **
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2019 on 26 degrees of freedom
Multiple R-squared: 0.4452, Adjusted R-squared: 0.3598
F-statistic: 5.216 on 4 and 26 DF, p-value: 0.003204
par(mfrow = c(2,2))
plot(lm)